% file to create graphs for Scruggs project
[labels,year,y1]=readColData('uk2.dat',10,0);
[labels,year,y2]=readColData('uk2tax.dat',10,0);
uibenef=y1(:,1);
uibeneftax=y2(:,1);
uidur=y1(:,2);
uiwait=y1(:,3);
is=y1(:,4);
istax=y2(:,4);
u=y1(:,5);
udur=y1(:,6);
waits=y1(:,7);
waitstax=y2(:,7);
atest0=y1(:,8);
atest1=y1(:,9);
%uidur2=y1(:,9);
%uibenef2=y1(:,10);
%uibenef2tax=y2(:,10);
tax=istax./is-1;
theta=[0.314 0.166 0.237 0.251 0.150 0.182 0.324 0.196 0.303 0.121 0.107 0.087 0.103 0.127 0.111 0.107 0.080 0.191 0.170 0.105 0.129 0.097 0.212 0.215 0.211 0.204 0.199 0.210 0.195 0.190 0.187 0.185]; 
thetatax=[0.391 0.225 0.313 0.325 0.229 0.286 0.426 0.271 0.392 0.199 0.185 0.117 0.133 0.187 0.241 0.247 0.110 0.171 0.150 0.135 0.169 0.137 0.142 0.135 0.171 0.174 0.149 0.160 0.135 0.130 0.067 0.105];
thetataxuidur=[0.391 0.225 0.313 0.325 0.229 0.286 0.426 0.271 0.392 0.199 0.185 0.117 0.133 0.187 0.241 0.247 0.110 0.171 0.150 0.135 0.169 0.137 0.142 0.135 0.171 0.164 0.099 0.100 0.135 0.070 0.067 0.095];
thetataxnoatest=[0.391 0.225 0.313 0.325 0.229 0.286 0.426 0.271 0.392 0.249 0.285 0.257 0.263 0.257 0.211 0.197 0.180 0.171 0.150 0.135 0.159 0.157 0.162 0.155 0.151 0.144 0.119 0.130 0.115 0.100 0.097 0.095];
thetanoatest=[nan nan 0.167 0.161 nan 0.202 0.144 0.146 0.143 0.141 0.307 nan nan nan nan 0.147 0.190 nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan]; %not corrected
thetataxtheta=[0.293 0.153 0.173 0.153 0.153 0.193 0.203 0.203 0.203 0.143 0.143 0.133 0.153 0.293 0.183 0.173 0.143 0.263 0.263 0.263 0.203 0.173 0.303 0.303 0.293 0.133 0.283 0.283 0.273 0.263 0.263 0.303]; 
thetataxunemp=[0.401 0.265 0.353 0.365 0.449 0.286 0.446 0.271 0.412 0.179 0.205 0.157 0.233 0.207 0.191 0.187 0.160 0.161 0.150 0.135 0.149 0.147 0.152 0.155 0.151 0.154 0.159 0.100 0.145 0.140 0.137 0.145];
thetataxunempdur=[0.491 0.255 0.453 0.225 0.459 0.286 0.506 0.261 0.442 0.199 0.185 0.107 0.263 0.167 0.231 0.227 0.220 0.151 0.140 0.135 0.139 0.127 0.172 0.175 0.121 0.104 0.099 0.100 0.145 0.090 0.087 0.095]; % not corrected OK now?
thetataxis=[0.401 0.235 0.313 0.325 0.219 0.266 0.446 0.251 0.402 0.189 0.155 0.097 0.113 0.257 0.141 0.247 0.230 0.171 0.150 0.135 0.149 0.157 0.162 0.155 0.151 0.144 0.119 0.130 0.115 0.100 0.097 0.095];
thetataxwdur=[0.401 0.225 0.313 0.315 0.229 0.286 0.426 0.271 0.392 0.189 0.175 0.117 0.133 0.187 0.241 0.247 0.110 0.171 0.150 0.135 0.169 0.137 0.142 0.135 0.171 0.174 0.149 0.160 0.135 0.130 0.067 0.105];
thetataxrho15=[0.391 0.225 0.283 0.275 0.219 0.286 0.406 0.461 0.422 0.219 0.215 0.147 0.163 0.197 0.191 0.187 0.200 0.161 0.140 0.125 0.139 0.147 0.152 0.145 0.141 0.134 0.109 0.120 0.105 0.090 0.077 0.085];
figure(1);
set(gcf,'DefaultLineLineWidth',1.5);
%plot(year,100*waits,'k-',year,100*(uidur.*uibenef+uidur2.*uibenef2)/(uidur+uidur2),'k--',year,100*is,'k-.',year,100*atest,'k:');
plot(year,100*waits,'k-',year,100*uidur,'k--',year,100*is,'k-.',year,100*(atest0+atest1)/2,'k:');
ylabel('%');
grid on;
title('Program benefit parameters, as % of average income');
legend('immediate benefits','UI benefits','income support','asset test');
print -depsc progparam;
figure(2);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*waitstax,'k-',year,100*uibeneftax,'k--',year,100*istax,'k-.',year,100*(atest0+atest1)/2,'k:');
ylabel('%');
grid on;
title('Program benefit parameters, as % of average income, with tax benefits');
legend('immediate benefits','UI benefits','income support','asset test');
print -depsc progparamtax;
figure(3);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,uiwait,'k-',year,uidur,'k--');
ylabel('weeks');
grid on;
title('Program eligibility parameters, in weeks');
legend('waiting period','UI eligibility period');
print -depsc progweek;
figure(4);
set(gcf,'DefaultLineLineWidth',1.5);
%plot(year,100*theta,'k-',year,100*thetatax,'k--');
plot(year,100*thetatax,'k-',year,100*uibenef,'k--');
ylabel('%');
grid on;
title('Effective UI coverage in simulated economies');
%legend('without tax benefits','with tax benefits');
legend('simulated benefits','UI benefits');
print -depsc figtheta;

%calculation of naive measures
beta=0.999165^0.5;
years=length(year);
naivetheta=zeros(years,1);
naivethetatax=zeros(years,1);
for i=1:years
  for j=1:2*uiwait(i)
    naivetheta(i)=naivetheta(i)+beta^(j-1)*waits(i);
    naivethetatax(i)=naivethetatax(i)+beta^(j-1)*waitstax(i);
    end;
for j=2*uiwait(i)+1:2*(uiwait(i)+uidur(i))
  naivetheta(i)=naivetheta(i)+beta^(j-1)*uibenef(i);
    naivethetatax(i)=naivethetatax(i)+beta^(j-1)*uibeneftax(i);
    end;
for j=2*(uiwait(i)+uidur(i))+1:100000
  naivetheta(i)=naivetheta(i)+beta^(j-1)*is(i);
    naivethetatax(i)=naivethetatax(i)+beta^(j-1)*istax(i);
    end;
%for j=2*(uiwait(i)+uidur(i))+1:2*(uiwait(i)+uidur(i)+uidur2(i))
%  naivetheta(i)=naivetheta(i)+beta^(j-1)*uibenef2(i);
%    naivethetatax(i)=naivethetatax(i)+beta^(j-1)*uibenef2tax(i);
%    end;
%for j=2*(uiwait(i)+uidur(i)+uidur2(i)+1):100000
%  naivetheta(i)=naivetheta(i)+beta^(j-1)*is(i);
%    naivethetatax(i)=naivethetatax(i)+beta^(j-1)*istax(i);
%    end;
naivetheta(i)=naivetheta(i)*(1-beta);
naivethetatax(i)=naivethetatax(i)*(1-beta);
  end;
figure(5);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*naivetheta,'k-',year,100*theta,'k--');
ylabel('%');
grid on;
title('Comparison of simulated and naive effective UI coverage');
legend('naive','simulated');
print -depsc fignaivetheta;
 
figure(6);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*naivethetatax,'k-',year,100*thetatax,'k--');
ylabel('%');
grid on;
title('Comparison of simulated and naive effective UI coverage with tax benefits');
legend('naive','simulated');
print -depsc fignaivethetatax;

figure(7)
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*u,'k-',year,udur,'k--');
ylabel('%, weeks');
grid on;
title('Labor market parameters');
legend('unemployment rate','unemployment duration');
print -depsc figlabormarket;

figure(8)
  set(gcf,'DefaultLineLineWidth',1.5);
xlaffer=0:0.001:1;
ylaffer=xlaffer-xlaffer.*xlaffer;
xlabel('\theta');
ylabel('value');
plot(xlaffer,ylaffer,'k-',xlaffer,0.24,'r-');
print -depsc figlaffer.eps;

figure(9);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*thetataxnoatest,'k-',year,100*thetatax,'k--');
ylabel('%');
grid on;
title('Robustness tests on UI coverage on asset test');
legend('no asset test','benchmark');
print -depsc figrobustnoatest;

figure(10);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*thetataxunemp,'k-',year,100*thetatax,'k--');
ylabel('%');
grid on;
title('Robustness tests on UI coverage: unemployment rate set to average value');
legend('counterfactual','benchmark');
print -depsc figrobustunemp;

figure(11);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*thetataxtheta,'k-',year,100*thetatax,'k--');
ylabel('%');
grid on;
title('Robustness tests on UI coverage: UI benefits set to average value');
legend('counterfactual','benchmark');
print -depsc figrobusttheta;

figure(12);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*thetataxunempdur,'k-',year,100*thetatax,'k--');
ylabel('%');
grid on;
title('Robustness tests on UI coverage: unemployment duration set to average value');
legend('counterfactual','benchmark');
print -depsc figrobustunempdur;

figure(13);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*thetataxwdur,'k-',year,100*thetatax,'k--');
ylabel('%');
grid on;
title('Robustness tests on UI coverage: wait time set to average value');
legend('counterfactual','benchmark');
print -depsc figrobustwdur;

figure(14);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*thetataxuidur,'k-',year,100*thetatax,'k--');
ylabel('%');
grid on;
title('Robustness tests on UI coverage: benefit duration set to initial value');
legend('counterfactual','benchmark');
print -depsc figrobustuidur;

figure(15);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*thetataxis,'k-',year,100*thetatax,'k--');
ylabel('%');
grid on;
title('Robustness tests on UI coverage: income support set to average value');
legend('counterfactual','benchmark');
print -depsc figrobustis;

figure(16);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*thetataxrho15,'k-',year,100*thetatax,'k--');
ylabel('%');
grid on;
title('Robustness tests on risk aversion');
legend('\rho=1.5','\rho=2.5 (benchmark)');
print -depsc figrobustrho;

figure(17);
set(gcf,'DefaultLineLineWidth',1.5);
plot(year,100*theta,'k-',year,100*thetatax,'k--');
ylabel('%');
grid on;
title('Robustness tests on tax benefits');
legend('no benefits','with benefits (benchmark)');
print -depsc figrobusttax;


%x=[ones(32,1) waits (uibenef.*uidur+uibenef2.*uidur2)./(uidur+uidur2) is atest uiwait uidur+uidur2 u udur];
%xtax=[ones(32,1) waitstax (uibeneftax.*uidur+uibenef2tax.*uidur2)./(uidur+uidur2) istax atest uiwait uidur+uidur2 u udur];
x=[ones(32,1) waits uibenef is (atest0+atest1)/2 uiwait uidur u udur];
xtax=[ones(32,1) waitstax uibeneftax istax (atest0+atest1)/2 uiwait uidur u udur tax];
[b,bint]=regress(theta',x);%'
impact=squeeze(max(x)-min(x))'.*b;%'
disp([b bint impact]);
[b,bint]=regress(thetatax',xtax);%'
impact=squeeze(max(xtax)-min(xtax))'.*b;%'
disp([b bint impact]);
%s=regstats(theta',x,'linear', {'beta','tstat','fstat','rsquare'});
%disp([s.beta s.tstat]);
%disp('R2, F');
%disp([s.rsquare s.fstat]);
